function [K j] = findgp (xy, grd)
 
   ## usage: [K j] = findgp (xy, grd)
   ## 
   ## find closest gridpoint

   global NETCDF

   if ischar(grd)

      if strcmp(NETCDF, "octcdf") && false # currently not safe
	 pkg load octcdf
	 nc = netcdf(grd, "r") ;
	 x = ncname(ncdim(nc){2}) ;
	 y = ncname(ncdim(nc){3}) ;
	 xx = nc{x}(:) ; yy = nc{y}(:) ;
      else
	 pkg load netcdf
	 for j=1:length(ncinfo(grd).Variables)
	    if (Vl=length(ncinfo(grd).Variables(j).Size)) > 2, break ; endif
	 endfor
	 VN = VJ = cell(1,Vl) ;
	 [VN{:}] = ncinfo(grd).Variables(j).Dimensions.Name ;
	 [VJ{:}] = ncinfo(grd).Variables(j).Dimensions.Unlimited ;
	 VJ = cellfun(@(v) v == false, VJ) ;
	 x = VN{1} ; y = VN{2} ;
	 xx = ncread(grd, x) ; yy = ncread(grd, y) ;	 
      endif
      [yy xx] = meshgrid(yy(:), xx(:)) ;
      if strcmp(x, "lon")
	 k = findgp(xy, [xx(:)'; yy(:)']) ;
      else
	 k = findgp(xy, [yy(:)'; xx(:)']) ;
      endif
      N = size(xx) ;
      [K(1,:) K(2,:)] = ind2sub(N, k) ;
      return ;

   else

      for j=1:columns(xy)
	 d = sumsq(grd - xy(:,j)) ;
	 K(j) = find(d == min(d))(1) ;
      endfor

   endif


endfunction
